#### no lead donors only -- uses only symmetric W matrix

rm(list=ls())
library(epicalc)
zap()


library(foreign)
library(MASS)

sar.f <- function(es, x, y, w){
    rho <- es[1]
    s2 <- es[2]^2
    if (s2<0) return(-1e20)
    k <- ncol(x)
    b <- as.matrix(es[3:(k+2)])
    llik <- 0
    nd <- nrow(w)
    n <- nrow(x)
    d <- n/nd
    y <- as.matrix(y)
    
    for(i in 1:d){
        A <- diag(1, nd, nd)-rho*w[,,i]
        dA <- det(A)
        if (dA <=0){
            llik <-  -1e20 
            break
        } else {
            aux <- log(dA) -nd/2 * log(s2)- (1/(2*s2)) * t(A %*% y[((i-1)*nd+1):((i-1)*nd+nd),] - as.matrix(x[((i-1)*nd+1):((i-1)*nd+nd),])%*%b) %*% (A %*% y[((i-1)*nd+1):((i-1)*nd+nd),]- as.matrix(x[((i-1)*nd+1):((i-1)*nd+nd),])%*%b)
            llik <- llik + aux
        }
    }
    return(llik)
}

sar <- function(x, y, w, stval,vars=c(""), routit=100000){
    mod <- optim(stval, sar.f, hessian = TRUE, method = "BFGS", 
            control = list(fnscale = -1, REPORT = 10, 
            reltol = 1e-16, trace = 1, 
            maxit = routit), x=x, y=y, w=w)
    vacov <- as.matrix(solve(-1*mod$hessian, tol=1e-24))
    se <- as.matrix(sqrt(diag(vacov)))
    zstat <- as.matrix(mod$par / se)
    pvalues <- round(2*pnorm(abs(zstat),lower.tail=FALSE),4)
    
    cat("results in order: rho, s, b \n")
    cat("parameter estimates: ", mod$par, "\n")
    cat("standard error: ", se, "\n")
    cat("z-statistics: ", zstat, "\n")
    cat("p-values: ", pvalues, "\n")
    list(llik=mod$value, par=mod$par, se=se, zstat=zstat, 
        pvalues = pvalues, 
        hessian=mod$hessian, vacov=vacov
        )
}

#routines for Clarke test (panel-wise)
sar.clr <- function(pars, x, y, w){
    rho <- pars[1]
    s2 <- pars[2]^2
    k <- ncol(x)
    b <- as.matrix(pars[3:(k+2)])
    nd <- nrow(w)
    n <- nrow(x)
    d <- n/nd
    y <- as.matrix(y)
    llik <- rep(NA, d)
for(i in 1:d){
    A <- diag(1, nd, nd)-rho*w[,,i]
    dA <- det(A)
    llik[i] <- log(dA) -n/2 * log(s2)- (1/(2*s2)) * t(A %*% y[((i-1)*nd+1):((i-1)*nd+nd),] - x[((i-1)*nd+1):((i-1)*nd+nd),]%*%b) %*% (A %*% y[((i-1)*nd+1):((i-1)*nd+nd),]- x[((i-1)*nd+1):((i-1)*nd+nd),]%*%b)    
    }
    return(llik)
}


setwd("C:/Users/Martin Steinwand/Documents/Project Lead Donor/Data")

rawdat <- read.dta("DonorRecipientDyads.4.2.dta")
#ensure data is sorted on r_ccode year d_ccode:
if (rawdat$r_ccode[1]!=rawdat$r_ccode[2]) cat("\n\n WARNING, data not sorted properly!!!\n\n")
if (rawdat$d_ccode[1]==rawdat$d_ccode[2]) cat("\n\n WARNING, data not sorted properly!!!\n\n")
if (rawdat$year[1]!=rawdat$year[2]) cat("\n\n WARNING, data not sorted properly!!!\n\n")

########################
#No lead donors only
workdat <- data.frame(rawdat[!rawdat$d_ccode==732 
                        & rawdat$ldstrong==0                        
                        & is.na(rawdat$lODAdis_f)==0 & is.na(rawdat$population)==0
                        & is.na(rawdat$voteshare_f)==0 & is.na(rawdat$rgdpch)==0
                        & rawdat$year >=1980 & rawdat$year <=1989,] )
attach(workdat)

vars <- cbind(AllPop, 1, ODADonYear*10e-2, hhi, OilExp_rf*1e-8, TotImp_rf*1e-8, 
        Camerica,Samerica, Africa, Meast, Asia, Oceania, population*1e-8, 
        rgdpch*1e-2, colown, voteshare_f)
    
vnames <- c("Aid", "Constant", "Total Donor Aid ", "Donor Concentration", 
                                "Oil Exports", "Total Imports", 
                                "Central America", "South America", "Africa", 
                                "Middle East", "Asia","Oceania", "Population", 
                                "GDP per capita, ppp","Former Colony", 
                                "Joint UN Voting")
colnames(vars) <- vnames

#get aggregates for 
nvars <- aggregate(vars, list(d_ccode, r_ccode), mean)
colnames(nvars)[1:2] <- c("d_ccode","r_ccode")
colnames(nvars)[ncol(nvars)] <- c("Joint UN Voting")
y <- nvars[,3]
x <- nvars[,4:ncol(nvars)]

#counts 
n <- nrow(nvars)
nd <- length(table(nvars$d_ccode))
d <- n/nd


#W matrix -- even
#array
wall.ev <- array(NA, c(nd, nd, d))
wraw <- matrix(c(rep(1/(nd-1),nd)),nd, nd, byrow=TRUE)
diag(wraw) <- 0
wall.ev[,,1:d] <- wraw


set.seed(021175)
stv <- runif((ncol(x) + 2))*.1
pop.no.wev <- sar(x, y, wall.ev, stval=stv, vars=vnames[2:length(vnames)])

setwd("C:/Users/Martin Steinwand/Documents/Project Lead Donor/Analysis/R&R.IO")
save(pop.no.wev, file="pop.pooled.1980s.wev")

#check model convergence, estimate with other starting values
#set.seed(072275)
#stv <- runif((ncol(x) + 2))*.1
#avg.no.wev2 <- sar(x, y, wall.ev, stval=stv, vars=vnames[2:length(vnames)])
